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ABSTRACT 

The Orion BN/KL complex is the nearest site of ongoing high-mass star formation. 
Recent proper motion observations provide convincing evidence of a recent (about 
500 years ago) dynamical interaction between two massive young stellar objects in 
the region resulting in high velocities: the BN object and radio Source I. At the same 
time, Source I is surrounded by a nearly edge-on disc with radius ^ 50 au. These 
two observations taken together are puzzling: a dynamical encounter between multi- 
ple stars naturally yields the proper motions, but the survival of a disc is challenging 
to explain. In this paper we take the first steps to numerically explore the preferred 
dynamical scenario of Goddi et al., in which Source I is a binary that underwent a 
scattering encounter with BN, in order to determine if a pre-existing disc can sur- 
vive this encounter in some form. Treating only gravitational forces, we are able to 
thoroughly and efficiently cover a large range of encounter parameters. We find that 
disc material can indeed survive a three-body scattering event if 1) the encounter is 
close, i.e. BN's closest approach to Source I is comparable to Source Fs semi-major 
axis; and 2) the interplay of the three stars is of a short duration. Furthermore, we 
are able to constrain the initial conditions that can broadly produce the orientation of 
the present-day system's disc relative to its velocity vector. To first order we can thus 
confirm the plausibility of the scattering scenario of Goddi et al., and we have signif- 
icantly constrained the parameters and narrowed the focus of future, more complex 
and expensive attempts to computationally model the complicated BN/KL region. 

Key words: binaries :close-circumstellar matter~ISM:individual objects (Orion 
BN/KL)-methods:A''-body simulations-stars:formation-stellar dynamics 



1 INTRODUCTION 

The details of how massive stars form are poorly understood. 
High-mass young stellar objects (YSOs) evolve rapidly and 
are intrinsically rare (due to the declining initial mass func- 
tion), which generally implies large distances (> 500 pc). 
Furthermore, they form deeply embedded in crowded re- 
gions where high extinction prevents observations in the 
optical and infrared. Lacking firm observational support, 
models of high-mass star formation have remained contro- 
versial, with three different scenarios active in the recent 
literature: (1) monolithic collapse and disc acc retion (e.g. 
lYorke fc SonnhalteJ [200^ : iMcKee fc Tan' '2003"): (2) com- 
petitive accretion in protoclusters (e.g. BonncU et al. 2001, 
[2OOJ), and similar theories in which protostellar dynam - 
ics greatly affect the ongoing accretion (iPeters et al.ll2010l) : 
and (3) stellar collisions and mergers jMurrav fc Linlll996l : 
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iBonnell et "allll998l : iBallv fc Zinneckeill2005l '). although this 
path is likel y only viable in the most extreme clustered en- 
viron ments (|Moeckel fc Clark3l201ll : lBaumgardt fc KlessenI 

As originally formulated, the last two models are based 
on the evidence that most massive s tars form in dense proto- 
stellar clusters (|Lada fc LadallioOSl l. where the stellar den- 
sities may be so high that dynamical interactions (includ- 
ing close passages or collisions and mergers) are naturally 
expected and may be common in the very early stages of 
massive cluster evolution. In contrast, the first scenario envi- 
sions a scaled-up version of low-mass star formation, where 
discs and outflows are invoked to reduce the effect of ra- 
diation pressure from massive protostars and enable ongo- 
ing accretion; discs concentrate the infalling material into 
small solid angles and coUimate the radiation field prefer- 
entially in the polar direction, where stellar photons can 
escape al ong the cavity of lower d e nsity gas excavated b y 
outflows (|Yorke fc Sonnhalteij[2003 : lKrumholz et al.l 120091 '). 
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Figure 1. Proper motions of BN and Source I in tlie ONC rest- 
fram e measured from m ulti-epoch radio continuum observations 
from [Coddi et aL ] l l201ll'l . The arrows indicate proper motion di- 
rection and displacement for 200 yr. The dash-dot line indicates 
motion backward over 560 years, assuming it is linear. The dashed 
lines indicate la uncertainties in position angles. The minimum 
(past) separation on the sky between Source I and BN was 50±100 
AU. The image in the inset shows the nearly edge-on compact 
disc around Source I as observed in 2006 with the VLA at 7mm 
ijGoddi et al.ll201lh . We note that another object in the region, 
source n, has been proposed to take pa rt in the same dyna mical 
interaction along with BN and Source I l| Gomez et al .1120081) . The 
analysis he re does not conside r source n, based on the evidence 
reported in lCoddi et al.l l|201lh . 



In recent years some modest convergence between model- 
ers has occurred, and the debate has advanced slightly. For 
instance, few today would argue that discs play no role in 
mediating accretion onto a massive YSO, or that radiation 
is unimportant. Whether the matter that ends up in a mas- 
sive star is bound prior to collapse as a massive core or is 
subsequently gathered from its environment remains a point 
of vigorous discussion, however, as does the importance of 
(proto) stellar dynamics. In order to constrain different the- 
ories, two different aspects should be investigated: on scales 
of individual YSOs, the properties of disc-outflow interfaces, 
in particular the structure of discs and the acceleration and 
coUimation mechanisms of outflows; and on scales of clus- 
ters, the role of dynamical interactions among protostars in 
the early phases of protocluster evolution. In these respects, 
the Orion BN/KL region can provide an ideal laboratory. 

The Orion BN/KL complex contains the nearest region 
of on going high-mass star formation (~420 pc jMenten et all 
l200'iD and is a compelling target for studying how massive 
YSOs form and interact with their surroundings. A dense 
protostellar cluster lies within the region containing two ra- 
dio sources that are believ ed to be massive YSOs: the highly 
embe dded radio Source I (iReid et al.ll2007l: [Matthews et al.l 
I2OI0I) and the BN object l|Becklin fc Neugebauej IiQGtI v^ 



which is the bright est source in the region in the mid- 
infrared at 12.4 /xm (jCezari et al.|[l998l ). Source I is proba- 
bly the best case known of a massive protostar with ongo- 
ing disc-mediated accretion. High angular resolution mon- 
itoring of molecular masers enabled for the first time the 
resolution of the launch/coUimation region of an outfiow 
from a compact disc, as well as tracing the dynamical evo- 
luti on of molecular g as over time at radii of 10-100 au 
(.Matthews et al.ll2010l ). The disc mass is poorly constrained, 
but it is est imated to lie in th e range 0.002 Mq< Mdisk < 
0.2 Mgfsee lGoddi et al.ll20l"ll '). Complicating matters, how- 
ever, recent works report strong evidence that a dynamical 
interactio n occurred 500 ye a rs ago between So urce I and 
BN (e.g. iGomez et all I2OO8I : ICoddi et al.]|201ll ). The stel- 
lar interaction probably resulted in high stellar proper mo- 
tions for Source I and BN (12 and 26 km s~^, respectively; 
iGoddi et allboill ) and possibly prod uced a fast bullet out- 
flow traced by CO and H2 emission (|Ballv et al.ll201ll ). In 
figure [l]we show the proper motions of Source I and BN, as 
w ell as the disc aroun d Source I; for complete details, refer 
to lGoddi et al] (|201lD . 

The strong evidence of a dynamical interaction and 
the simultaneous existence of a circumstellar disc (radius 
~ 50 AU) around Source I are challenging to explain. Two 
possibilities have been considered: an existing disc is de- 
stroyed in the encounter and then rebuilt in the following 
500 years from material dispersed from the original discs 
or interstellar gas in the region, via Bondi-Hoyle accretion 
(jBallv et al.l I2OIII ): alternatively, the disc is truncated at 
the radius of the close passage but survives the encounter 
(jGoddi et al.ll20l"ll ). The former authors speculate about the 
magnetoyhydrodynamic consequences of a rapidly changing 
potential due to chaotic stellar dynamics. The latter au- 
thors focused on the stellar dynamics of the stars alone, 
performing n-body experiments to explore different dynam- 
ical scenarios, settling on a preferred model of an interaction 
between an existing binary and a single star. 

While the full story of the complicated BN/KL region 
remains unclear, steps forward must involve modeling of the 
system. Any modeling of the system's dynamics must be per- 
formed from a set of initial conditions that permits a stellar 
interaction to: (1) result in high post-encounter velocities of 
the interacting stars, and (2) preserve some of the original 
circumstellar material despite chaotic and potentially com- 
plicated stellar dynamics. We turn to numerical experiments 
to attempt to constrain this set. In this work, we take a first- 
order approach to the problem, considering only the most 
important physical process: the gravitational interaction be- 
tween the stars and its effect on the circumbinary material. 
This is reasonable, as the system's mass is dominated by the 
stars while the gas around Source I is most likely rotationally 
supported. A fully realistic approach should include, beyond 
gravitational dynamics, hydrodynamic, radiative, and mag- 
netic forces. In that case, however, parameter space would 
be tremendously expanded, and likewise the computational 
cost of a single experiment. The final explanation of BN/KL 
will undoubtedly involve an interplay between all of these 
forces; we are not presenting this work as a full explanation 
of the region, but rather as the first step in a systematic, 
step-by-step approach to the problem. 

In this paper, for the first time we numerically investi- 
gate effect of a gravitational scattering encounter between 
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Figure 2. Schematic representation of an encounter of a single 
star with a binary in the binary's rest frame. The gray curved line 
shows the path of the hyperbolic encounter the binary and the in- 
truder would experience in the two-body approximation where the 
binary is represented by its center of mass. We show the param- 
eters we sample over (black labeled quantities), as well as some 
auxiliary variables that are a consequence of the initial conditions 
(gray labeled quantities). The periastron separation between the 
binary center of mass and the intruder rp, their relative velocity 
at infinity v^o, and the impact parameter b are related by the 
gravitational focusing relationship shown in equation [2] The an- 
gles 9 and <j> determine the angle between the angular momenta 
of the binary-single orbit, Li,_s, and the orbit of the binary L^j^, 
and p determines the initial orbital phase of the binary. 



stars on a circumbinary disc. The paper is structured as fol- 
lows. The details on simulations are described in ^ Main 
results from the numerical work are presented in !j3l In iJH we 
discuss these results in the context of existing observations 
of Orion BN/KL, and in Sj5]we summarize our conclusions 
and look forward to future steps. 



2 SIMULATION DETAILS 

In this work we consider an isolated n-body system com- 
prising of three stars, arranged as an initial binary and 
a single star, and a circumbinary disc, which we model 
with low-mass particles. Hydrodynamic forces in the disc 
are neglected, an approximation which allows us to com- 
plete many thousand scattering experiments in a reason- 
able timeframe. Large numbers of experiments are neces- 
sary, due to the chaos inherent to the three-body problem. 
It is impossible to definitively trace backwards from the cur- 
rent system, observed with uncertainty, to an earlier con- 
figuration; instead statistical studies must be performed to 
determine the likelihood of similar outcomes to what we ob- 
serve. TlM2i^od2_system is evolved using Aarseth's nbody6 
code l|Aarsethll2000l . l2003l '). which utilises a fourth-order Her- 
mite inte grator. While the code includes a KS regularization 
scheme (jKustaanheimo &: Stiefej Il965h for handling close 
encounters between bodies, because of the atypical setup 
of this problem (three dominant bodies and a disc of low- 
mass particles), we suppress this feature and the equations 
of motion are integrated in standard Cartesian space. Be- 
cause NBODY6 does not smooth the force between particles, 
there is in principle the possibility that undesired relaxation 



effects can occur in the disc; however in practice, over the 
short timescale of these simulations this is not a concern. 

The two main elements of this problem are the in- 
teraction between the disc and the stars, and the grav- 
itational scattering encounter between the stars. Sep- 
arately, each of these elements have been extensively 
explore d: the t heory of binary-single scattering is laid 
out in Heggid (1 19751) . and a long history of nu m erical 
work i ncludes [Seckerl (ll92ol). ISaslaw et al.l (1 19741 ) iHilld 
(1 19751). IValtonen fc Heggi3 ~ (| 19791 ). iHut fc Bahcalll (| 19831 ). 
and iHud (|l993l )." Disruptive disc-star encounters have like- 
wise been studied numerically in a va r iety of contexts 
(e.g. IClark c fc Prinaie 'l993'; HcUcr 19951; iHall et al.l 1 19961 : 
iBoflin et al..,1998. : .Moeckel fc Bally, 20oi ). The combination 
of the two processes is unstudied to our knowledg^Ui a-nd 
interesting in general. Parameter space is vast, and it is for- 
tunate that BN/KL provides a physical reference from which 
to take initial steps. 



2.1 Initial conditions 

We ran two sets of simulations, which we label in terms of 
the masses of the initial stellar configuration in Mq. One 
set of 4500 simulations, (10,3) + 10, has a low-mass member 
in the binary, while the other set of 7000 runs, (10,10)-|-10, 
has all three stars with equal masses. Each experiment be- 
gins with a circular binary with semi-major axis a = 10 au 
surrounded by a circumbinary disc with total mass 0.1 Mq. 
The results of this paper are, strictly speaking, only applica- 
ble to discs that are in this non-self gravitating, dynamically 
unimportant regime. While the dynamics of the actual scat- 
tering interaction between the stars is probably insensitive 
to the disc mass, a more massive disc could potentially alter 
the trajectory of the stars in the immediate leadup to the 
stellar interactions. Note that the estimated remnant disc 
mass of less than 0.2 Mq is comfortably unimportant to 
the dynamics of a massive binary over short timescales. The 
simulated disc is trunc ated at the inner edge at the 3:1 reso- 
nance with the binary IjArtvmowicz fc Lubowll 19941 ) which is 
approximately 2.08a = 20.8 au, extends to the outer edge at 
10a = 100 au, and follows an arbitrary but plausible surface 
density profile E(r) oc . Because we ignore hydrodynamic 
forces there is no need to resolve the disc vertically, and the 
disc particle are initially coplanar with the binary. The disc 
is composed of 2048 particles, and tests with 1024 and 4096 
particles show that the results are essentially independent 
of disc resolution. 

In order to place the (initially unbound) third star, the 
binary is treated as a single body at its center of mass, and 
a Keplerian orbit is determined for some relative velocity 
at infinity Voa and desired periastron Vp. This periastron is 
the unperturbed separation that occurs when the binary is 
treated as its center of mass; the actual orbits of the three 
bodies diverges from this path as they approach each other 
and the two-body approximation to the system breaks down, 
but this serves to define the initial conditions. The third star 
is introduced at a separation of 110a = 1100 au from the 



^ Although there has been some work in a similar spirit, inves- 
tigating the surviva l of accretion discs in black hole interactions 
llLin fc Safilawll 19771) . 
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Figure 3. Examples of a successful interaction {top row) and an unsuccessful one {bottom row) from the (10,10)+10 runs. Left column: 
the paths of the three stars through the interaction in a frame drifting with respect to the center of mass of the system. The initially- 
single star is the darker line. Numbers with arrows label the stars at the beginning and end of the encounter. Right columns: the final 
distribution of the stars (black points) and disc particles (gray points) 500 years after the scattering event. Disc particles considered part 
of the disc by the criteria discussed in the text are plotted as larger blue points. We show both a wider view in the center of mass frame, 
and closer view centered on the final binary. 



binary. At this separation, the ratio of the perturber's force 
to the binary's force at the disc edge is less than 10~^. 

In the two-body approximation the orbit is determined 
by just Voo and Vp, but there are also several angles re- 
quired to identify a set of initial conditions, specifying the 
orientation of the binary. These can be described by the 
angle between the angular momentum of the binary's orbit 
(T-ibin) and the angular momentum of the binary-single or- 
bit (Lf,_s), given by the polar angle 6 and azimuthal angle 
(j}; the orientation of the binary's Laplace-Runge-Lenz vector 
relative to some reference vector in the binary plane; and the 
mean anomaly of the binary, M. The Laplace-Runge-Lenz 
vector is a conserved vector that describes the orientation 
of the binary's ellipse, and the mean anomaly is related to 
the position in the orbit. Because we consider only circular 
binaries in this work, these are degenerate and can be re- 
placed by the orbital phase p of the binary, measured from 
the reference vector in the binary plane. For the reference 
vector we use Lb_s x hbin- In figure [2] we show a schematic 
representation of our initial conditions 

Because our primary goal is not to dete rmine the prob- 
abilit y of a certain scattering event (as in e.g. iHut fc Bahcalll 
1 19831 ) but rather to see what encounter parameters can re- 
sult in disc retention, we have sampled uniformly over the 
periastron of the initial orbit rather than the impact pa- 
rameter squared. This means that the initial conditions we 
choose are not formally proportional to their probability, 
but it ensures that we have good coverage over potentially 
interesting regions of the Vp-Voo plane that could otherwise 
be sampled with low probability by a standard scattering 



experiment setup. However, as we show below, the two sam- 
pling methods are virtually equivalent over the parameter 
range we cover. 

2.2 Parameter ranges 

Our initial conditions for each experiment are then described 
by five random numbers sampled with a quasi-random num- 
ber generator, namely: 

rp/au e [0, 60) 
Woo /km s"' e [0, 15) 
cos e G [-1,1) 
■^e [0,27r) 
PG [0,27r). 

We have given the limits on Vp and Voo in physical units, but 
they can also be expressed in non-dimensional form relative 
to the natural length and velocity scales of the problem. 
Reasonable choices for these are the binary semi-major axis 
a, yielding rp/a G [0, 6), and the critical velocity of the three- 
star system, Vc- This is the value of Voo at which the three- 
body energy is zero; if mi and m2 are the masses of the 
binary components and ma is the intruder, it is given by 

G mim2(mi + m2 + ma)] ^''^ 

Vc ~ —, T ■ (1) 

_ a m3(mi + m2) J 

For the (10,10)-|-10 runs we have Vc — 36.5 km s~^ and 
Voo/vc G [0,0.41), and for the (10,3) + 10 runs we have — 
21.7 km s"^ and Uoo/wc G [0,0.69). 
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The upper limit for r-p is determined by what might be 
an 'interesting' encounter, that is the relative velocity be- 
tween the binary and the single is boosted during the three- 
body interaction. This requires a close triple approach and 
hardening of the binary; previous work suggests that the 
maximum interaction distance between the bi nary and the 
single is a few times t he binary separation fe.g. lSaslaw et aU 
I1974I : IValtone n"l974Y Our limit of 60 au, or 6a, safely cov- 
ered the full range of potentially interesting periastra. After 
an initial run of 2000 experiments, we determined that no 
interesting encounters were taking place beyond — 30 
au for the (10,3) + 10 case. After this an additional 1500 
experiments were performed, drawn from the same quasi- 
random sampling series but carrying out only those runs 
with rp ^ 30 au. The (10,10)-|-10 case had a few interesting 
encounters at larger periastra, so for this set all 5000 ex- 
periments were performed over the full range. The limit on 
Voo is set by velocity dispersion of the Ori on Nebula Clus- 
ter, a ~2-3 km s"^ l|jones fc Walkeij|l988l '). If the velocity 
dispersion in the region is Maxwellian, the integrated prob- 
ability of an encounter with Vao > 15 km s~^ is less than 1 
per cent with a ID velocity dispersion of 3 km s~^. 

As noted above, strictly speaking the initial conditions 
are not chosen according to equal probabilities. Such a sam- 
pling would be proportional to the area of the surface ele- 
ment orthogonal to Voo, and is achieved by a uniform sam- 
pling of where b is the impact parameter. However, the 
underlying physical motivation for this work is such that 
the stellar masses, range of interesting periastra, and likely 
encounter velocities place us in a gravitationally focused en- 
counter regime, where the paths of the stars deviate sub- 
stantially from rectilinea r motion (see figure [ 2]). Fo llowing 
the familiar derivation in lBinnev fc Tremaind ( 20081 ) . in the 
two-body approximation the periastron separation between 
the binary and the single is related to the masses, impact 
parameter and velocity at infinity via 



b = rr, 



2G(mi +m2 + ma) 



(2) 



The second term in brackets is associated with the gravita- 
tional focusing of the orbits. When this term is dominant, 
the right hand side is approximately linear in Vp, in which 
case b^ « rp to a good approximation. The largest rp we use 
is 60 au, and the smallest value of the focusing term occurs 
with Voo = 15 km s~^ in the (10,3)-|-10 runs, where it takes 
a value of approximately 181 au. Thus our entire parameter 
range is in a strongly focused regime, and our uniform sam- 
pling of rp is very close to the equal-probability sampling, 
which is formally proportional to b^. We can use our setup 
as a reasonable proxy for one specifically tailored to deter- 
mining scattering outcome cross sections, and we do this in 
Appendixl^to verify that our results are in agreement with 
previous work. 



2.3 Experiment termination and characterisation 

The goal of these experiments is to determine if there are 
regions in parameter space that can lead to the retention of 
disc material in a binary-single encounter that has a large 
final relative velocity. In general, an encounter between a 
binary and a single star can result in one of two outcomes: 
ionisation of the system into three single stars, or an end 



state consisting of a binary and a single. The latter possibil- 
ity can be further divided into a flyby, in which the initially 
single star remains single, or an exchange in which one of 
the original binary members is displaced by the intruder. For 
our choice of initial conditions ionisation is energetically for- 
bidden, as the internal energy of the binary is greater than 
the kinetic energy associated with the binary-single orbit, 
i.e. Voo < Vc- All of our experiments must terminate in a 
binary and a single. 

The interplay of three bodies in a scattering event can 
be quite complex a nd last for an unpredictable length of 
time (e.g. figure 3of lH^ll993f) . and some way to determine 
the end of a simulation is required beyond integrating for a 
fixed period. During each experiment we track the size of the 
three stars' configuration, measured by their total moment 
of inertia in the center of momentum frame, 



We terminate the runs when two criteria are met. First, in 
the two body approximation (i.e. when the binary is treated 
as its center of mass) the binary-single orbit is hyperbolic. 
Second, 500 years have past since the last minimum in I3, 
which we then define as the time of the scattering. At this 
point we have the three stars arranged as a binary and a 
single, each component surrounded (or not) by a cloud of 
disc particles. In physical units, the integrations range from 
about 10^ to 10* years, roughly following a power law as 
between those limits. 

In order to determine final remnant disc membership, 
we rely on a circularisation sc heme simil a r to one that has 
been used in previous work (|Halj 1 19971 : iMoeckel fc Ballvl 
I2OO6I ). At the end of the simulation, the binary is treated 
as a single body at its center of mass. Disc particles that 
are energetically bound to this center of mass have passed 
the first cut. These then have their orbital angular momen- 
tum calculated, and a Keplerian orbital period determined. 
In other applications, the disc particles may be placed on 
circular orbits corresponding to their specific angular mo- 
mentum, the argument being that viscous effects and shocks 
between flows of material on crossing orbits will damp any 
eccentricity. In this way a circularised and settled disc can 
be reconstructed without simulating hydrodynamic forces 
for long time periods. 

Because of the specific constraints of the present prob- 
lem, some care must be taken and we do not directly re- 
construct the disc profile in this way. Instead, we impose 
two additional cuts to determine the final disc membership. 
First, the orbital period of the disc particle about the binary 
must be less than 500 years, the time between the scatter- 
ing event and the present day. Second, the periastron of the 
particle's orbit around the binary center of mass must be 
less than 50 au, the observed size of the disc. These dual 
constraints ensure that the disc particle will have passed 
within 50 au at least once during the 500 years available 
to circularise the remnant disc, and will probably have had 
several orbital periods to settle into a disc-like state. We do 
not attempt to reconstruct the disc profile at all, but rather 
just track the amount of disc material that satisfies these 
constraints and thus estimate the mass of the disc within 50 
au. 

This method removes any tidal tails or wispy extended 
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clouds from our calculation of the remnant disc mass. In 
some cases it excludes material that could very well be part 
of a the disc, either on longer timescales or in a full hydrody- 
namic calculation, but because of the tight time constraints 
of this problem we do not assign it to the remnant. Note 
that current observations of the disc (of both SiO masers 
and radio continuum) show no signs of extended structure 
at radii larger than 50 au, although future sensitive obser- 
vations of thermal molecular lines and dust continuum (e.g. 
with ALMA) may reveal weaker emission associated with 
extended gas remnants. 

Our approach should be conservative; including hydro- 
dynamics, so that particles are not free to pass through the 
disc but would rather shock and become incorporated into 
the disc, would probably increase the amount of material in 
the disc. With the mass of the disc remnant interior to 50 au 
determined, we deem an encounter successful if it results in 
both a high velocity v^ei between the binary and the single, 
and retains at least 10 per cent of the initial disc around 
the binary. In figure [S] we show examples of two encoun- 
ters with equal mass stars, one simple (and successful) and 
one complex, showing the paths of the stars and the final 
distribution of disc particles 500 years after the encounter. 
The eccentric appearance of the disc in the simple case is at 
least partly a consequence of the dissipationless physics we 
are modeling; the detailed morphology of the disc remnants 
should not be interpreted as meaningful in this work. 



2.4 Integration accuracy 

The distribution of energy errors in the simulations is ap- 
proximately log-normal, with the mean log(Ai5/ii'o) = -4.93 
and standard deviation 0.53. Experiments with the binary 
and its circumbinary disc in isolation show that the accu- 
mulated error in the integration of the disc particles and 
the stellar particles are of the same order of magnitude. If 
we were concerned about the precise details of a particu- 
lar set of initial conditions, we would strive for higher ac- 
curacy. Since we are more interested in a broad statistical 
survey of parameters, previous work has argued that more 
lenient energy tolerances are acceptable as long as conclu- 
sions d on't rely on the details of a few integrations. I ValtonenI 
found little accuracy dependence in the statistics of 
a suite of approximately 200 simulations with mean rel- 
ative ene rgy errors ranging f rom 5 x 10~* to 3 x lO"'^. 
IPflamm-A ltcnburg & Kroupal (|2006l ) compared the statisti- 
cal results of 1000 four-body decay experiments with mean 
energy errors ranging from about 6 x 10~ to 6 x 10~ , 
and recovered the same results to within lo-. lHut fc Bahcalll 
l|l983ll used a relative energy error li mit of 10~^, and thei r 
results compare very well to those of iFregeau" et all (|2004 ). 
In Appendixl^we show that our experiments, although not 
tailo r ed to the same goal, acceptably match the results of 
[Hu3 (|l993l ). We conclude that our energy tolerance is suffi- 
cient. 



3 RESULTS 

3.1 Initial conditions resulting in success 

We begin by identifying which regions of parameter space 
can lead to a successful scattering event. Recall that we de- 
fine success as an experiment in which a relative velocity Vrei 
between the single and the binary of at least 20 km s~^ is 
obtained, while at the same time 10 per cent of the original 
disc material is retained around the final binary. 



3.1.1 Results for (10, 3)+10 

This set of 4500 encounters begins with an unequal mass bi- 
nary composed of a 10 Mq primary and a 3 Mq secondary. 
In the left side of figure |4] we show four projections through 
the space of initial conditions, covering the most interest- 
ing combinations of parameters. In each panel the graylevel 
shows the percentage of runs in each bin that were success- 
ful. 

Looking first at the Vao-Tp plane, we note that all suc- 
cessful encounters have a periastron less than roughly 20 au. 
Because of this, after completing 2000 runs we restricted 
the remaining 2500 to have r-p < 30 au. Initial conditions 
were drawn from the same quasi-random number genera- 
tor as all the other runs, we simply did not integrate them 
if they missed our cut. The sampling density over the 
interesting portion of parameter space is thus identical to 
the (10,10)-|-10 set, consisting of 7000 experiments. When 
constructing the histograms that follow, the small number 
of points with large velocities at rp > 30 au are weighted 
according to their lower sampling density, although this is 
a minor effect. There were some \ow-Voo encounters with 
larger that resulted in high velocities, although none of 
them retained a disc. These high-r^ encounters all occured 
near cos 9= 1, a nearly co-planar prograde encounter. If the 
phase of the binary orbit is favorable, this configuration can 
maximize the interaction time between the intruder and the 
secondary. The panels on the right show that there is no 
preferred combination of Voo and cos 6 leading to success- 
ful encounters, and no preferred orientation of the intruder's 
initial orbit. 

In the left side of figure [S] we show the marginal dis- 
tributions of the encounter results for some of the parame- 
ters (i.e. all other parameters have been integrated over). In 
three of the panels we plot histograms showing the fraction 
/ of encounters with high final v^i ('fast') and successful 
encounters ('fast -f disc') for bins of Vp, Voo, and cos 0. The 
fast-|-disc histograms have Poisson la error bars. There is 
a clear tendency for low Voo cases resulting in more fast 
results, but this results in only a mild dependence for the 
fast-|-disc runs, concentrated in the lowest Voo bin. There is 
likewise structure in the cos 6 histogram for the fast cases, 
but the fast-|-disc histogram is virtually featureless. In con- 
trast, there is an obvious dependence on rp, with a dramatic 
falloff in successful encounters over the range < Vp < 25. 
Small-r-p encounters are not only more likely to result in a 
high ejection velocity, but the fraction of those fast encoun- 
ters that retain a disc likewise increases with smaller Vp. 
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Figure 4. Four projections through parameter space for the (10,3)+10 case (left panels) and the (10,10)+10 runs (right panels). The 
graylevel shows the percentage of initial conditions resulting in both a high final relative velocity between the binary and the single, 
v,.^l > 20 km s~^, and retention of 10 per cent of the disc around the binary. In both sets, small values of the periastron separation rp 
are preferred for successful encounters. The equal mass runs have a higher percentage of successful encounters with prograde (cos 9 > 0) 
encounter orientations. 
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Figure 5. Histograms for the (10,3)+10 runs (left panels) and (10,10)+10 runs (right panels) showing the fraction of runs that resulted 
in iVe! > 20 ("fast"), and the fraction that also retained 10 per cent of the initial disc ("fast+disc") as a function of rp, Voa, cos 9, and 
Nmin(^3). The last is the number of minima in the three-body moment of inertia 1^ showing the distributions of encounter complexity 
for those same sets of runs. 



3.1.2 Results for (10, 10)+10 

In this set of 7000 encounters all three stars have the same 
mass, 10 Mq. We again begin by looking at four projections 
through the initial conditions, shown on the right side of 
figure[4l In contrast to the (10,3) + 10 case there are encoun- 



ters at large periastron resulting in success, those occurring 
at low and prograde inclinations. Besides the increased 
range of rp that allows disc retention and the somewhat 
higher percentages at low rp, the most striking difference be- 
tween these runs and the case with a 3 Mq secondary is the 
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fraction of disc remaining 

Figure 6. The percent of runs with high velocities (A^^) that also 
retain a disc {Nf^_^) for different fractional disc retention criteria. 
The results quoted in the bulk of the paper are for 10 per cent 
disc retention. 



angular dependence of successful encounters. The density of 
successful initial conditions increases significantly with pro- 
grade encounters (0 <cos 6< 1) relative to retrograde cases. 

On the right side of figure[5]we show the marginal distri- 
bution histograms for the equal mass encounters. The gen- 
eral trend of the rp dependence is the same as the (10,3)4-10 
case, although with a broader range. While there is some 
bias toward low-Hoo cases resulting in high final velocities, 
this trend is again less apparent in the fast-|-disc cases. The 
cos 6 histogram contrasts sharply to the (10,3)-|-10 experi- 
ments. At retrograde inclinations the fraction of successful 
runs is fairly flat, and begins to rise as the encounters be- 
come prograde. 



3.1.3 Dependence on the remnant disc mass cutoff 

Our choice of 10 per cent of the disc remaining bound as a 
criterion for success is somewhat arbitrary. Given the lack of 
constraints on both the initial and current disc masses, this 
is in some sense unavoidable. This choice does have some 
effect on the results. In figure [§] we show the fraction of 
encounters with Vrei > 20 km s~^ that retain discs of dif- 
ferent masses. As the cutoff on the remaining disc mass in- 
creases fewer runs result in success, with a precipitous drop- 
off above about 20 per cent of the disc remaining, especially 
for the equal-mass runs. 

To a good approximation, the behavior of the results 
shown in figures [4] and [5] scale with these curves. Versions 
of these figures with cutoffs of 2.5 and 20 per cent rather 
than 10 per cent display very little functional variation in 
the initial conditions leading to successful outcomes, and a 
simple renormalisation suffices for the purposes of this pa- 
per. Thus our quoted results with the 10 per cent criterion 
can be thought of as representative of results up to a re- 
quirement of ~ 20 per cent disc retention, modulo a scaling 
of order unity; above this, the order of magnitude of the 
results changes and success becomes much less likely. 



(10,3)+10 (10,10)+10 




20 30 40 50 60 70 80 90 20 30 40 50 60 70 80 90 

Vrel / km S"^ V,e| / km s"' 

Figure 7. The distribution of final velocities v^^i for the 
(10,3)-f 10 case (left) and the (10,10)-f 10 case (right). Shown are 
the runs that resulted in v^ei > 20 km ("fast"), the subset 
that also retained 10% of the initial disc ("fast+disc"), and for 
the (10,3)-|-10 runs the subset that retained the lightest star as 
a member of the final binary ("fast-|-disc-|-3"). For the (10,3)+10 
plot, all cases with v^^i > 85 km s~^ are represented in the final 
bin. 

3.2 Disc survival with encounter complexity 

Intuitively, one would expect a complex encounter (such 
as that shown in the lower row of figure ^ to disrupt a 
disc much more efficiently than a relatively clean and fast 
exchange. This intuition is borne out by the experimental 
results. In the lower-right panels of figure [S] we show the 
number of close triple encounters, measured by the number 
of minima Nmin in I3, for all the experiments performed, 
for the high-Vrei runs, and for successful encounters. These 
plots are normalized to the total number of experiments per- 
formed for each set. The distribution of N,nin is quite broad, 
peaked at low values but extending from 10s to 100s of min- 
ima. 

The fraction of all runs that are classified as fast is fairly 
independent of Nmin, at around 10 per cent of all encoun- 
ters. In contrast, the encounters that leave a remnant disc 
are almost all single-minimum events. For the (10,3)-|-10 case 
less than 3 per cent of successful encounters go through mul- 
tiple close triple passages, and for (10,10)4-10 runs less than 
6 per cent do. This is the most restricting component of a 
successful scattering event; it must be a clean encounter, ei- 
ther exchanging the intruder for a binary member or ejecting 
the intruder with a minimal amount of interaction. 

3.3 Final relative velocities and binary orbits 

In figure [7] we show the distribution of the final relative ve- 
locities Vrei for both sets of runs. There are several previous 
results concerning the escape velocity of the decay of three 
equal-mass bodies from an initi ally binary-free setup, e.g . 
iMikkola fc ValtonenI (|l986l ) and ISterzik fc DurisenI (|l99a ). 
with good agreement to an analytic approximation gener- 
alized from iHeggid (Il975l ). These distributions rise sharply 
to a maximum and decay over a longer range of velocities. 
While not directly applicable to our study, the general be- 
havior of our equal-mass results mirrors the high-velocity 
tail of those earlier works. The (10,3)4-10 results show more 
structure than the equal-mass case. Here we plot both the 
distribution of all encounters that retain a disc, as well as 
those that keep the 3 Mq star in the binary. These cases are 
exclusively below v^ei = 35 km s~^, with all the higher ve- 
locity cases resulting from scatterings where the 3 Mq star 
is ejected. At higher relative velocities the fast and fast-|-disc 
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Figure 8. The location of the final binaries in the semi-major axis— eccentricity plane for all cases with v^ei > 20 km s~^. Graylevels 
show the density of points in each cell. The histograms show the distribution integrated along the eccentricity, i.e. showing the fraction 
of runs in bins of a only. For the (10,3)-|-10 runs, the cases that retain the 3 Mq star are shown in the black histrogram, while those 
where it is ejected are shown in grey. 



populations are nearly identical. This is due to momentum 
conservation in the three-body scattering. These encounters 
are all ones in which the 3 Mq star is swapped for a 10 
Mq star and the recoil velocity of the now more massive bi- 
nary is much less than that of the liberated low-mass star. A 
low binary velocity combined with a deeper potential makes 
disc retention easier. 

The velocity behavior is reflected in the orbit of the bi- 
nary after the scattering, which we show in figure [S] The 
relative velocity between the binary and the single at the 
end of the scattering experiment is a function of the change 
in energy of the binary as a result of the encounter, and any 
increase in Vrei compared to the initial Voo must be compen- 
sated by an increase in the binary's binding energy. For the 
equal mass cases considered in the (10, 10) -I- 10 runs the en- 
ergy change is simply inversely proportional to the binary's 
semi-major axis; those runs with Vrei > 20 km s~^ leave 
binaries with a semi-major axis of 6.7 ± 1.0 au and eccen- 
tricity 0.66 ± 0.22. For simplicity we have characterised the 
distributions by their means and standard deviations, al- 
though the eccentricity distributions in particular are quite 
broad. While these equal mass stars are effectively identi- 
cal in these simulations, we note that in 51 percent of the 
encounters with «,ei > 20 km s~^ an exchange occurs, and 
the initially single star takes the place of one of the original 
binary members. 

With the unequal masses in the (10,3)4-10 runs, the fi- 
nal binary's semi-major axis depends on which star is ejected 
from the system. When the 3 Mq star is retained in the bi- 
nary, the semi-major axis must be reduced in order to boost 
the relative velocity; the semi-major axis distribution for 
these runs is characterized by a mean value of 5.2 ± 0.77 



au, with eccentricities 0.78 ± 0.18. When the 3 Mq star 
is ejected, however, the semi-major axis of the now more- 
massive binary can be larger than the initial value. The 
semi-major axis distribution for these runs is much broader 
than the other cases. There is a main peak around 12.6 ±5.6 
au and eccentricities 0.67 ± 0.24, but there is also a popula- 
tion above 20 au with high eccentricity. This population is 
problematic for the survival of a roughly co-planar remnant 
disc of 50 au. 

The high eccentricity of all the encounters leading to 
high velocities is notable; if Source I is found to be a cir- 
cular binary, it seems unlikely that it took part in a recent 
scattering event. Note that the binarity of Source I is not 
observationally established, as length scales < 10 au are be- 
yond the resolution of current instruments. 



3.4 Remnant disc orientation 

A final observable in the BN/KL system that we can com- 
pare our simulations to is the relative orientation of the 
disc and the presumed binary's velocity vector. The proper 
motion of Source I is from north-west to south-east, and 
the nearly e dge-on disc has a spin axis aligned north-east- 
south- west (|Goddi et al.ll201ll . see figure [T]). We can search 
our data for experiments that yield this roughly orthogonal 
alignment. In figure [9] we show for both sets of simulations 
the relationship between two angles for all successful runs: 
the angle between the disc and the binary angular momen- 
tum, cos~^(Ldiac • Lbin), and the angle between the disc and 
the binary-single relative velocity, cos~^(Lciisc • Vrei). The 
disc angular momentum is taken to be the sum of the an- 
gular momenta of all disc particles that are deemed part of 
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Figure 9. The large panel of each set of plots shows the distribution of successful encounters in the plane of two disc orientation angles; 
the angle between the remnant disc orientation and the binary's orbit, and the angle between the disc orientation and the binary's 
velocity. For these purposes the disc and binary orientations are taken as their respective angular momentum vectors. Observations have 
shown that the disc symmetry axis and the proper motion of Source I are roughly orthogonal (see Figure 1). We have overlaid boxes 
around the region of this plane where these vectors are within 20 degrees of orthogonality. The blue box encompasses results where the 
binary and disc angular momenta are within 30 degrees of alignment, and the red box encompasses all other binary— disc alignments. 
The smaller plots trace these results back to their input parameters, showing the density of initial conditions in the rp-cos 8 plane that 
lead to final disc orientations in the red and blue regions of interest. 



the final disc remnant, calculated in the reference frame of 
the binary. While each set has some higher-density regions 
in this plane, in general most combinations of these angles 
are represented. 

We have highlighted all runs that have a mostly orthog- 
onal orientation between the disc angular momentum and 
Vrel, which we take as within 20 degrees of orthogonal, with 
boxes. The blue box encompasses all nearly-orthogonal runs 
that also have the disc and the binary aligned to within 30 
degrees; the red box encompasses all other disc-binary ori- 
entations. In the smaller panels of the figure, we trace back 
to the initial conditions to see if these interesting boxes are 
populated from a specific region of parameter space, plot- 
ting the density of runs that end up populating each box. 
The only interesting features turns out to be in the r^-cos 
6 plane. The first notable result is that the (10,3)-f 10 series 
contains no runs that led to alignment between the binary 
and a disc that is orthogonal to the relative velocity be- 
tween the single and the binary. The equal mass series, in 
contrast, has several runs in this category. Crucially, the ini- 
tial conditions these runs came from are much smaller than 
the total parameter range we covered. Roughly 10 per cent 
of all initial conditions with periastra less than about 20 au 
and prograde encounter angles result in high stellar veloc- 
ities and remnant discs with similar orientations to Source 
I. 



3.5 Likely effect of relaxing assumptions 

All of our simulations have taken the binary as initially cir- 
cular. We can draw on extensive binary-single scattering 
work to es timate the effect of this simplifying assumption. 
[Hu^ (|l993l ') presents scattering results for equal-mass sys- 
tems at three eccentricities: e = 0, 0.7, and 0.99. He shows 
that the cross sections for scattering results (see the ap- 
pendix for some discussion of this) depend only weakly on 



the eccentricity; increasing the eccentricity from to 0.7 
approximately doubles the cross section for both simple ex- 
change encounters and resonant encounters, and a further 
increase to 0.99 has minimal effect. This increase in en- 
counter cross section is mainly manifested in enc ounters tha t 
only slightly change the binary energy, however (|Hutlll984l '). 
Over the range of velocities we consider, the cross section for 
a given energy change to the binary as a result of a scatter- 
ing encounter in the range of interest, i.e. those encounters 
that generate Vrei > 20 km s~^, are almost identical for all 
three eccentricities. The total number of sufficiently velocity- 
boosting encounters is thus largely insensitive to the eccen- 
tricity of the binary, and our zero-eccentricity assumption 
should not affect these results. 

We also assume that the binary and its disc are largely 
dynamically unprocessed, in the sense that the disc is copla- 
nar to the binary. Our results do rely on this assumption; if 
the disc and the binary are seriously misaligned prior to the 
encounter, our effort to trace back the interesting outcomes 
to their initial conditions (figure [9| is compromised. While 
there would almost certainly be a set of encounters that are 
more likely to lead to the proper disc-«,.e; orientation, it 
might not be as nicely contained as in the coplanar case. 

The initial parameters of our disc, for which we have at- 
tempted to choose plausible values, are also an assumption. 
If the original disc is considerably larger or more compact 
than our choice of 100 au, the fractional amount of disc ma- 
terial retained will be different than the values quoted here. 
Because the bulk of the retained material is originally in the 
inner regions of the disc, the conclusion that non-negligible 
amounts of material can remain bound should hold true as 
long as the inner regions of the disc are still intact. The 
effect of tidal or direct interactions with low mass cluster 
members prior to the encounter with the BN object should 
not greatly alter these results either, provided they do not 
get captured into a wide binary /hierarchical triple; although 
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in this case the disc w ould eventually be fully disrupted 
iMoeckel fc Ballvll2006l ). and this seems incompatible with 
the observed system. 

Finally, we are working under the assumption that the 
gravitational interaction between the stars and the disc is 
the dominant physical mechanism at play. We are in effect 
simulating the dense midplane of the disc, rather than the 
surface or atmosphere that is actually observed in radio con- 
tinuum and SiO masers. Our assumption is thus that ma- 
jority of the disc mass is in the molecular gas of the disc 
midplane, obscured from observations. If the disc were to- 
tally ionised, ignoring the hydrodynamic forces would cer- 
tainly not be a valid approximation. Other effects that may 
work to sculpt the disk morphology, such as photoevapo- 
ra tion, could be at work around a hot source like Source 
I l|Hollenbach et al]|l994l ). Indeed, with the nominal rem- 
nant disc mass in our simulations of ~ 0.01 Mq, a pho- 
toevaporative wind with M ~ 10~^ Mq yr~^ could act 
rapidly enough to be important over 500 yr. We have worked 
mostly with the fractional disc mass rather than its absolute 
value, however; increasing our initial disc mass by an order 
of magnitude, which would still be quite small compared to 
the binary mass, would make the orbital timescale the only 
relevant one over 500 yr. While we regard the 'gravity as 
the dominant force' approximation to be a reasonable one 
for shaping the size of the remnant disc, further work with 
more massive discs (including at least hydrodynamics) is 
necessary to confirm this. 



4 DISCUSSION 

Disc survival in encounters involving more than two stars 
is an unexplored problem. The combination of star-disc 
interactions and three body scattering is quite relevant 
in star formation theories where stellar dynamics play a 
dominant (or at least perturbing) role in the ac cretion 
process fe.g. iBonneU et aLlliooH : IPeters et al.l I2OI0I '). It is 
notable that both Orion BN/KL and Cepheus A HW2 
ICunningham et al.ll200sj ). two of the closest and best stud- 
ied sites of ongoing massive star formation, show evidence 
of stellar interactions during the accretional growth of a 
massive YSO. The parameter space of binary-single scat- 
tering with discs present is huge, and the strong evidence 
of a multiple encounter in the BN/KL region has provided 
a physically motivated point to focus our initial investiga- 
tions around. Our simulations described in Section [3] enable 
us to constrain the likely encounter parameters that could 
lead to the observed BN and Source I velocities and disc 
configuration. 

The first point we make is that a scenario with all three 
bodies of similar mass seems more likely than one with a 
low mass object. The low-mass star is more likely to be 
the ejected member of a scattering encounter; in our exper- 
iments, about 80 percent of encounters resulting in high ve- 
locities end up with the two most massive stars in the binary. 
If these mass ratios are similar to reality, this would imply 
that the BN object is likely the lowest mass star in the trio, 
which is incompatible with luminosity-bas e d mass estimates 
of ar ound 8-12 Mq l|Scoville et al.ril983l : iRodn'guez et al.l 
I2OO5I ). No similar estimate of Source I's mass ex ists, as it is 
undetected in the mid-IR tareenhill et"alll2004l ). However, 



based on the motion of the SiO masers and the ionising flux 
derived from radio continuum emission, it is est imated to 
be greater than 7-10 Mq fsee iGoddi et al.ll201ll . for a fufl 
discussion of Source I's mass). It is possible that the scat- 
tering encounter was one of the 20 percent that leaves the 3 
Mq star in the binary. In this case, the overall magnitude of 
the relative velocity between Source I and BN could be con- 
sistent with the observed value, and the final binary would 
be guaranteed to have a smaller semi-major axis after the en- 
counter than before. However, momentum conservation and 
the observed proper motions suggest a mass ratio of Source 
I to BN of approximately a factor of two (assuming that the 
center of mass of the three stars is stationary in the Orion 
rest frame). The remnant disc is a further line of evidence 
arguing against a low mass member of the system. While 
retaining disc material is not any more of an issue for the 
(10,3)4-10 cases than the equal mass runs, the orientation of 
the disc relative to the proper motion vector is problematic. 
Not a single run ended up generating the observed angles, 
as discussed in section 13.41 We thus turn our attention to 
the equal mass case. 

Referring to figure [5] the total percentage of encounters 
that result in both high final velocities and a disc is not large 
(<10 percent). Importantly, though, ^50 percent of those 
('fast') encounters also retain 10 percent of the initial disc 
('fast-|-disc'). Interestingly, the distribution of successful en- 
counters is heavily biased toward small encounter distances; 
almost 30 percent of encounters with rp < 5 au result in 
success. This is somewhat at odds with the intuitive picture 
many would have, that the current disc size marks the clos- 
est approach of the intruder, as this truncation mechanism 
has been noted and investig ated in many numerical studies 
of st a r-disc encounter s (e.g. Clarke fc Pringld 19931: Hellei 
I995I : iHall et al.l Il996l : iBoflBn et al.lll998l : iMoeckel fc Bally 
20061 ). If the high proper motions were not present, trunca- 



tion would in fact be an attractive explanation for the disc 
size. However, as seen in figure O encounters well within the 
current size of the disc are necessary to achieve binary hard- 
ening and the consequent relative velocity boost. In fact, no 
encounters outside 50 au, the observed size of the disc, result 
in velocity boosts. 

The key to reconciling very close encounters and disc 
survival lies in the distribution of encounter complexities. 
All successful encounters are over very quickly, mostly with 
a single minimum in the spatial configuration of the bodies. 
The picture to have in mind is a clean interaction taking 
place within the circumbinary cavity, essentially a head on 
encounter between the binary and the intruder, that sends 
the stars on their new paths with minimal direct disruption 
of the disc. Thus rather than the surviving disc being a 
truncated remnant after the binary and the intruder undergo 
a Keplerian encounter, as in previous work with star-disc 
encounters, it consists of the material that is dragged along 
as the binary's trajectory is almost impulsively altered. 

This scenario has several advantages compared to the 
idea of accreting the disc from the ISM after the encounter. 
While there is some evidence that small short-lived discs 
can be set up around an accreting object moving through 
material wit h some momentum gradient perpendicular to 
its velocity l|Ruffertlll997l , Il999l ). the disclike structures in 
those simulations are smaller as a fraction of the gravita- 
tional radius (r^ < O.lra) than the one around Source I 
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(with re = 2GM/v^ ~ 250 au, compared to the 50 au 
disc). Furthermore, the time between the scattering event 
and now is only ~500 years and the distance traveled is 
only ~700 AU, i.e. a few times the gravitational radius. Set- 
ting up an identifiable Bondi-Hoyle type fiow over that short 
of a time and distance may be difficult. In the picture pre- 
sented here, the material starts out bound and in a disc, and 
immediately after the encounter the remnant is still largely 
co-planar and bound to the binary. Over the ~ 500 years 
since the encounter, this material only has to reorganize it- 
self from the eccentric orbits we see in our simulations into 
the disc observed today rather than build itself from scratch 
from the ISM. 

In summary, our numerical work has enabled us to con- 
strain the parameters of a stellar scattering that probably 
occurred in BN/KL: roughly equal mass stars, a prograde 
encounter between the binary orbit and the orbit of the 
binary-single system, and a periastron comparable to the bi- 
nary semi-major axis. We note that this work includes only 
the most dominant force, gravity. This choice provides the 
most computationally efficient strategy for large parameter 
searches. We cannot, however, say anything about the cir- 
cularisation of the remnant disc without including at least 
hydrodynamics. Additionally, magnetic fields ma y play a sig- 
nifica nt role in generating the observed outflows I Bally et al.l 
I2OIII ). as well as providing some support to the dis c around 
Sour ce I (see the discussion of this point in Goddi et al.l 
I2OIII ). Despite these limitations, this work has significantly 
narrowed the range of parameters for future studies which 
will incorporate more complicated and computationally ex- 
pensive physics. 



5 CONCLUSIONS AND FUTURE WORK 

We have shown that the gravitational dynamics of a velocity- 
boosting three-body scattering event between an existing 
binary with a circumbinary disc and a single star do not 
preclude the survival of some fraction of the disc. Despite 
the chaos of the three-body encounter and the nearly im- 
pulsive velocity kick the binary center of mass experiences, 
there are regions of encounter parameter space that permit 
at least 10 percent of the disc to remain associated with 
the binary. The encounter type that is most permissive of 
disc survival is a nearly head-on one, so that the stellar dy- 
namics occurs mainly in the cavity of the circumbinary disc. 
Encounters that are complex with a long-lived stellar inter- 
action will more efficiently disrupt the disc. When we take 
into account the present day orientation of the disc around 
Source I, we are able to further constrain the encounter to 
one with roughly equal-mass stars, and with a prograde ori- 
entation between the orbit of the binary and that of the 
binary-single orbit. 

We stress that this work is the first step in a system- 
atic investigation of the dynamics at work in BN/KL, hence 
the conclusions must be interpreted with the first-order ap- 
proach we took in mind. As we have noted throughout the 
paper, additional physics such as hydrodynamic, magnetic, 
and radiative forces will all be important to fully understand 
the BN/KL region. Looking forward, we cannot speculate 
too far about the path our explorations will take, as each 
addition of physics will build upon the results of the pre- 



vious work. The next step is clear, however, and that is to 
include hydrodynamics and study the settling of the disc in 
a self-consistent fashion, drawing initial conditions from the 
restricted set that we have identified here. This will address 
the greatest shortcoming of this study, which is the lack of 
dissipation in the disc. Given the 500 year time constraint 
of this problem, this is an important issue to investigate for 
this scenario, and is the subject of ongoing work using an 
n-body/SPH code. 
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of the binary-single system is zero (see equation [1} . For our 
setup, recall that Vc in physical units is 36.5 km s"'^. As our 
limits on Voo rule out ionisation of the binary, all of our simu- 
lations are in the regime Vao < Vc- W hile this rules o u t com - 
parison wit h the more recent tests in iFregeau et al] (|2004l ). 
figure 3a of |Hu3 (|l993l ) includes the case of an equal mass, 
zero eccentricity binary scattering with a star of matching 
mass over a velocity regime coincident with ours. 

As a benchmark we use the data from that figure, and 
for consistency with that paper we define an exchange as an 
encounter with fewer than four minima in I3 after which one 
member of the binary has been replaced by the intruder, and 
all encounters with more than four minima are deemed to 
be in resonance. In studies designed to study cross sections, 
rather than sampling the initial conditions evenly over Vp, 
one samples over the impact parameter b according to its 
probability, i.e. uniformly in b^. The total cross section for 
outcome X is estimated simply by 

/ \ I \ («oo ) , . , V 

ox [Vao ) = -nb^ax («oo ) -, T • ( Al) 

Ntot[Voa) 

Nx is the number of experiments resulting in the outcome 
of interest, and Ntot is the total number of experiments per- 
formed out to a maximum impact parameter bmax- While 
we have chosen to evenly sample rp instead, in the regime 
where Voo < Vc the gravitational focusing term of the en- 
counter cross section dominates, and rp oc fe^ to a good ap- 
proximation (see the discussion in section [23}. We thus use 
equation lAll directly to calculate cross sections, after bin- 
ning our data logarithmically in Voo- The errors associated 
with this t ype of Monte Carlo estimate are fully discussed in 
iHut fc Bahc all ( 1983) ; we show the statistical uncertainty 

Aax{v^) = crx{vac,)/\/Nxiv^)- (A2) 

Over most of our velocity range, we sample out to large 
enough Vp that no more exchanges or resonances are oc- 
curring, i.e. the encounters are all flybys. In this case we 
can regard our experiments as fully sampling the parameter 
space of interest. At the lowest velocities we may be miss- 
ing some very large b encounters that could contribute to 
the cross sections, but the Hut data ends before that point 
anyway and no comparison to previous work can be made 
below Voo/vc ~ 0.08. The comparison is shown in figure \Kl\ 
with the cross sections normalized to the area of the binary 
orbit with semi-major axis a. Although our error bars are 
slightly larger than the earlier work, which used about twice 
times as many encounters optimized to this calculation, the 
results of our experiments agree at about the la level. 

This paper has been typeset from a T^^X/ I^TfjX file prepared 
by the author. 



APPENDIX A: CONSISTENCY CHECK VIA 
CROSS SECTIONS 

As a point of comparison between our simulations and pre- 
vious work, we estimate the total cross section for exchange 
and resonance in the equal-mass case (10,10)-|-10. Cross sec- 
tions are customarily presented with Voo measured in terms 
of the critical velocity Vc, which is where the total energy 
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Figure Al. Cross section for direct excliange from our 
(10,1 0)+10 case (dots withi error bars), compared to the results 
from lHu3 l|l993l ) (connected lines) over the range in common with 
our parameters. The cross section for exchange is shown in black, 
and the cross section for a resonant interaction is shown in red. Er- 
ror bars from the Hut paper are not shown, but they are slightly 
smaller than ours. The gray dashed line included for reference 
and is the same one shown in the Hut paper, which is analytic 
expression for the high-velocity li mit of the i onisation cross sec- 
tion, given by o- = (40/9)7ra2/-y-2 | |Hutlll983l) . 
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